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We investigate the dynamics of relativistic spinning test particles in the spacetime of a rotating 
black hole using the Papapetrou equations. We use the method of Lyapunov exponents to determine 
whether the orbits exhibit sensitive dependence on initial conditions, a signature of chaos. In the 
case of maximally spinning equal-mass binaries (a limiting case that violates the test-particle approx- 
imation) we find unambiguous positive Lyapunov exponents that come in pairs ±A, a characteristic 
of Hamiltonian dynamical systems. We find no evidence for nonvanishing Lyapunov exponents for 
physically realistic spin parameters, which suggests that chaos may not manifest itself in the gravi- 
tational radiation of extreme mass-ratio binary black-hole inspirals (as detectable, for example, by 
LISA, the Laser Interferometer Space Antenna). 



PACS numbers: 04.70.Bw, 04.80.Nn, 95.10.Fh 



I. INTRODUCTION 

The presence of chaos (or lack thereof) in relativis- 
tic binary inspiral systems has received intense attention 
recently due to the implications for gravitational-wave 
detection HlllSHIElSIEIl, especially regarding the 
generation of theoretical templates for use in matched fil- 
ters. There is concern that the sensitive dependence to 
initial conditions that characterizes chaos may make the 
calculation of such templates difficult or impossible @- 
In particular, in the presence of chaos the number of tem- 
plates would increase exponentially with the number of 
wave cycles to be fitted. In addition to this important 
concern, the problem of chaos in general relativity has 
inherent interest, as the dynamical behavior of general 
relativistic systems is poorly understood. 

Several authors have reported the presence of chaos for 
systems of two point masses in which one or both par- 
ticles are spinning 0, 0, @ • Our work follows up on Q , 
which studies the dynamics of a spinning test particle or- 
biting a nonrotating (Schwarzschild) black hole using the 
Papapetrou equations [Eqs. l|2.7f) ]. We extend this work 
to a rotating (Kerr) black hole, motivated by the expecta- 
tion that many astrophysically relevant black holes have 
nonzero angular momentum. Furthermore, the potential 
for chaos may be greater in Kerr spacetime since the Kerr 
metric has less symmetry and hence fewer integrals of the 
motion than Schwarzschild. In addition, the decision to 
focus on test particles is motivated partially by the LISA 
gravitational wave detector Q, which will be sensitive 
to radiation from spinning compact objects orbiting su- 
permassive black holes in galactic nuclei. Using the Kerr 
metric is appropriate since such supermassive black holes 
will in general have nonzero spin. 

There are many techniques for investigating chaos in 
dynamical systems, but for the case at hand we favor the 
use of Lyapunov exponents to quantify chaos. Informally, 



if eo is the phase-space distance between two nearby ini- 
tial conditions in phase space, then for chaotic systems 
the separation grows exponentially [sensitive dependence 
on initial conditions): e(r) = eoe Ar , where A is the Lya- 
punov exponent. (See Sec. lIII Al for a discussion of issues 
related to the choice of metric used to determine the dis- 
tance in phase space.) The value of Lyapunov exponents 
lies not only in establishing chaos, but also in providing 
a characteristic timescale t\ = 1/ A for the exponential 
separation. 

By definition, chaotic orbits are bounded phase space 
flows with at least one nonzero Lyapunov exponent. 
There are additional technical requirements for chaos 
that rule out periodic or quasiperiodic orbits, equilibria, 
and other types of patterned behavior |Toj . For exam- 
ple, unstable circular orbits in Schwarzschild spacetime 
can have positive Lyapunov exponents but such or- 
bits are completely integrable (see Sec. I VI I and hence not 
chaotic. In practice, we restrict ourselves to generic or- 
bits, avoiding the specialized initial conditions that lead 
to positive Lyapunov exponents in the absence of chaos. 

The use of Lyapunov exponents is potentially dan- 
gerous in general relativity because of the freedom to 
re-define the time coordinate. Chaos can seemingly 
be removed by a coordinate transformation: simply let 
t' = logr and the chaos disappears. Fortunately, in our 
case there is a fixed background spacetime with a time 
coordinate that is not dynamical but rather is simply 
a re-parameterization of the proper time. As a result, 
we will not encounter this time coordinate re-definition 
ambiguity (which plagued, for example, attempts to es- 
tablish chaos in mixmaster cosmological models, until 
coordinate-invariant methods were developed [Tl|i. Fur- 
thermore, we can compare times in different coordinate 
systems using ratios: if t p is the period of a periodic orbit 
in some coordinate system with time coordinate t, and r p 
is the period in proper time, then their ratio provides a 
conversion factor between times in different coordinate 
systems 
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For chaotic orbits, which are not periodic, we use the 
average value of dt/dr over the orbit, so that 

as discussed in Sec. lVII Dl [This more general formula re- 
duces to Eq. in the case of periodic orbits.] Since we 
want to measure the local divergence of trajectories, the 
natural definition is to use the divergence in local Lorentz 
frames, which suggests that we use the proper time r as 
our time parameter. The Lyapunov timescale in any co- 
ordinate system can then be obtained using Eq. I|1.2|l . 

Lyapunov exponents provide a quantitative definition 
of chaos, but there are several common qualitative meth- 
ods as well, none of which we use in the present case, 
for reasons explained below. Perhaps the most common 
qualitative tool in the analysis of dynamical systems is 
the use of Poincare surfaces of section. Poincare sections 
reduce the phase space by one dimension by consider- 
ing the intersection of the phase space trajectory with 
some fixed surface, typically taken to be a plane. Plot- 
ting momentum vs. position for intersections of the tra- 
jectory with this surface thengives a qualitative view of 
the dynamics. As noted in Q, such sections are most 
useful when the number of degrees of freedom minus the 
number of constraints (including integrals of the motion) 
is not greater than two, since in this case the resulting 
points fall on a one-dimensional curve for non-chaotic or- 
bits, but are "dusty" for chaotic orbits (and in the case of 
dissipative dynamical systems lie on fractal attractors). 
Unfortunately, the system we consider has too many de- 
grees of freedom for Poincare sections to be useful. It is 
possible to plot momentum vs. position when the trajec- 
tory intersects a section that is a plane in physical space 
(say x = 0) 0|, but this is not in general a true Poincare 
section. 1 

Other qualitative methods include power spectra and 
chaotic attractors. The power spectra for regular orbits 
have a finite number of discrete frequencies, whereas their 
chaotic counterparts are continuous. Unfortunately, it is 
difficult to differentiate between complicated regular or- 
bits, quasiperiodic orbits, and chaotic orbits, so we have 
avoided their use. Chaotic attractors, which typically 
involve orbits asymptotically attracted to a fractal struc- 
ture, are powerful tools for exploring chaos, but their use 
is limited to dissipative systems |l(j| . Nondissipative sys- 
tems, including test particles in general relativity, do not 
possess attractors [l2 |. 



In |3, they are aided by the symmetry of Schwarzschild space- 
time, which guarantees that one component of the spin ten- 
sor (Sec. Ill Xl below - ) is zero in the equatorial plane. As a result, 
it turns out that all but two of their variables are determined 
on the surface, and thus their sections are valid. Unfortunately, 
the reduced symmetry of the Kerr metric makes this method 
unsuitable for the system we consider in this paper. 



Following Suzuki and Maeda Q, we use the Papa- 
petrou equations to model the dynamics of a spinning 
test particle in the absence of gravitational radiation. We 
extend their work in a Schwarzschild background by con- 
sidering orbits in Kerr spacctime, and we also improve 
on their methods for calculating Lyapunov exponents. 
The most significant improvement is the use of a rigor- 
ous method for determining Lyapunov exponents using 
the linearized equations of motion for each trajectory in 
phase space (Sec. 1111 A*fr . which requires knowledge of the 
Jacobian matrix for the Papapetrou system (Sec. IVB|l . 
We augment this method with an implementation of an 
informal deviation vector approach, which tracks the size 
of an initial deviation of size eo and uses the relation 
e(r) = e e Ar discussed above. We are careful in all cases 
to incorporate the constrained nature of the Papapetrou 
equations (Sec. Ill A|l in the calculation of Lyapunov ex- 
ponents (Sec. lIVBl . 

We use units where G = c = 1 and sign conventions as 
in MTW 0. We use vector arrows for 4- vectors (e.g., p 
for the 4-momentum) and boldface for Euclidean vectors 
(e.g., £ for a Euclidean tangent vector). The symbol log 
refers in all cases to the natural logarithm log e . 

II. SPINNING TEST PARTICLES 

A. Papapetrou equations 

The Papapetrou equations 01 describe the motion of a 
spinning test particle. Although Papapetrou first derived 
the equations of motion for such a particle, the formula- 
tion by Dixon 15] is the starting point for most inves- 
tigations because of its conceptual clarity. Dixon writes 
the equations of motion in terms of the 4-momentum p a 
and spin tensor S a @, which are defined by integrals of 
the particle's stress-energy tensor T a/3 over an arbitrary 
spacelike hypersurface E: 

p Q (£) = / T Q/3 cE/3 (2.1) 

S a/3 (z, £) = 2 / (x [Q - z^ a )T^~< d£ 7 , (2.2) 

where z is the coordinate of the center of mass. The 
equations of motion for a spinning test particle are then 



Vffp" = -^K a pVS^ (2.3) 
V^S"" = 2p^v v \ 

where v M is the 4- velocity, i.e., the tangent to the par- 
ticle's worldline. It is apparent that the 4-momentum 
deviates from geodesic motion due to a coupling of the 
spin to the Riemann curvature. 
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1. Spin supplementary conditions 



of motion are 



As written, the Papapetrou equations 12.31 are under- 
determined, and require a spin supplementary condition 
to determine the rest frame of the particle's center of 
mass. Following Dixon, we choose 



P^" = 0, 



(2.4) 



which picks out a unique worldline that we identify as the 
center of mass. In particular, in the zero 3-momentum 
frame defined by p 1 = 0, applying Eq. I|2.4jl to Eq. (|2.2jl 

yields 



T t x i T°° d 3 x 

J /.—const. 

L . Too d 3 x ■■ 

J t=const. 



(2.5) 



which is the proper relativistic generalization of the New- 
tonian center of mass. The frame defined by p l = is 
thus the rest frame of the center of mass, and in this 
frame Eq. (|2.4I) implies that S' -' = 0, i.e., the spin is 
purely spatial in the rest frame. 

A second possibility for the supplementary condition 



is 



v u S^ = 0. 



(2.6) 



This condition has the disadvantage that it is satisfied 
by a family of helical worldlines filling a cylinder with 
frame-dependent radius |l5l Il6| , centered on the world- 
line picked out by condition 12.41 As a result, we adopt 
PfiS^ = as the supplementary condition. 

We note that the difference between the conditions 12. 41 
and 12.61 is third order in the spin [which follows from 
Eq. (|2.13() below], which means that it is negligible for 
physically realistic spins (Sec. Ill B|) . In particular, the 
two conditions are equivalent for post-Newtonian expan- 
sions where condition ^. 61 is typically employed jig . 



2. A reformulation of the equations 

For numerical reasons, we use a form of the equations 
different from Eqs. I|2.3I) . (We discuss this and other 
numerical considerations in Sec. IV Ah Following the ap- 
pendix in [3} , we write the equations in terms of the mo- 
mentum 1-form p^ and the spin 1-form S^. 2 The system 
under consideration is a spinning particle of rest mass /1 
orbiting a central body of mass M; in what follows, we 
measure all times and lengths in terms of M, and we 
measure the momentum of the particle in terms of /i, so 
that p v p v = — 1. In these normalized units, the equations 



2 The lowered indices are motivated by the Hamiltonian formula- 
tion for a nonspinning test particle, where it is the one-form p M 



where 



Vv Pu = -Rl^v v p a S p (2.7) 

Vff Sfi = -Pfj, i.R* a p lS SaV^PjSs 



R*Y V = \R% a e^. (2. 



The tensor and vector formulations of the spin are re- 
lated by 



Su = 



ti — 2 



L V S aS3 



and 



Uf3, 



(2.9) 



(2.10) 



where u v = p v / 'fi (= p v in normalized units). In addition, 
the spin satisfies the condition 



C c/i _ 1 q qfiv _ 02 



(2.11) 



where S is the spin of the particle measured in units 
of \xM (see Sec.lTTRl). 

Because of the coupling of the spin to the Riemann 
curvature, the 4-momentum p^ [Eq. (|2.1|l ] is not paral- 
lel to the tangent v^. The supplementary condition 12.41 
allows for an explicit solution for the difference between 
them (see 0] for a derivation): 



where 



and 



W 



, = -* R *^ Sa ppS 1 



(2.12) 



(2.13) 



(2.14) 



The normalization constant N is fixed by the constraint 
v^v^ — — 1. We see from Eq. (|2.13|l that the difference 
between p M and is 0(S 2 ), so that the difference be- 
tween Eqs. andESJis 0(S 3 ). 

The spin 1-form satisfies two orthogonality constraints: 



and 



P^S U = 



= 0. 



(2.15) 



(2.16) 



that is canonically conjugate to x M 1 1 .'jj . 



These two constraints are equivalent as long as is 
given by Eq. (l2~l~2T) . since w^S^ oc * R*^^ S^S a = 0. 
When parameterizing the initial conditions, we enforce 
Eq. H2.15fl : since we use Eq. (|2.12|l in the equations of 
motion, Eq. I|2.16|l is then automatically satisfied. 
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3. Range of validity 

We note that the Papapetrou equations include effects 
due only to the mass monopole and spin dipole (the pole- 
dipole approximation). In particular, the tidal coupling, 
which is a mass quadrupole effect, is neglected. It is also 
important to note that the Papapetrou equations are con- 
servative and hence ignore the effects of gravitational ra- 
diation. For a thorough and accessible general discussion 
of the Papapetrou equations and related matters, includ- 
ing a comprehensive literature review, see Semerak |l9j. 



B. Comments on the spin parameter 

It is crucial to note that, in our normalized units, the 
spin parameter S is measured in terms of /iM, not /x 2 . 
The system we consider in this paper is a compact spin- 
ning body of mass fi orbiting a large body of mass M, 
which we take to be a supermassive Kerr black hole sat- 
isfying M k 1O 5 -1O 6 M . We will show that physically 
realistic values of the spin must satisfy S <C 1 for the 
compact objects (black holes, neutron stars, and white 
dwarfs) most relevant for the test particles described by 
the Papapetrou equations. 3 The case of a black hole is 
simplest: a maximally spinning black hole of mass [i has 
spin angular momentum s = /i 2 , so a small black hole /i 
orbiting a large black hole of mass M /i has a small 
spin parameter S: 



s 



< 



M 



< 1. 



The limit is similar for neutron stars: most models of 
neutron stars have a maximum spin of s max w 0.6 |20| . 
which gives S < 0.6 fi/M. 



1. Bounds on S for stellar objects 

The bound on S is relatively simple for black holes 
and neutron stars, but the situation is more complicated 
for compact stellar objects such as white dwarfs. The 
maximum spin of a stellar object is typically determined 
by the mass-shedding limit, i.e., the maximum spin be- 
fore the star begins to break up. The spin in the case 
of the break-up limit is the moment of inertia times the 
maximum (break-up) angular velocity: s max = If2 max . 
If we write I — a\xR 2 and r2 max = j3\J Gfi/R 3 for some 
constants a, (3 < 1, then we have 



= a{3{Gfi 3 R) 1/2 . 



(2-17) 



3 Recall that the Papapetrou equations ignore tidal coupling, so 
they are inappropriate for modeling more extended objects. 



The values of a and f3 depend on the stellar model; if 
we use the values for an n — 1.5 polytrope, we get 
a = 0.2044 and (3 = 0.5366 21], so that s max = 

o.no (<V^) 1/2 - 

The limit in Eq. (|2.17() depends on the mass-radius 
relation for the object in question. Since most neutron 
stars have masses and radii in a narrow range, the esti- 
mate of s max » 0.6 [i 2 discussed above is sufficient, but 
for white dwarfs the value of s max can depend strongly 
on the mass. An analytical approximation for the mass- 
radius relation for non-rotating white dwarfs is [22] : 4 



— = 0.01125 



, -1/3 
" * fM 1/2 



where 



and 



m = i 



4/3 



1.454 M, 



©• 



(2.18) 



(2.19) 



(2.20) 



We could plug Eq. (|2~T%1) into Eq. lt2~T7Jl to obtain an 
order-of-magnitude estimate, but [2l| tabulates a con- 
stant J equivalent to the product a(3 (which increases as 
the angular velocity of the star increases). They write 
J = J {GAP Rq) 1 / 2 for a rotating white dwarf, where 
J depends on the polytropic index n of a non-spinning 
white dwarf of the same mass, and i?o is the non-spinning 
radius. In our notation, this reads 



J(G/j 3 R) 



1/2 



(2.21) 



White dwarfs with // > 0.6 Mq are not well approximated 
by polytropes (the effective polytropic index varies from 
near 3 in the core to near 1.5 in the outer parts), but 
useful bounds can be obtained by substituting R from 
Eq. H2.18fl . which is more accurate for white dwarfs than 
a pure polytrope model. Plugging 12.181 into 12.211 and 
converting to geometric units gives 



= 77.68J// 3 A4 /3 /(M) 1/4 - 



(2.22) 



From Table 3 in |2l|, we have J — 0.1660 for a maxi- 
mally rotating n — 1.5 polytrope (vs. a/3 — 0.110 for a 
slowly rotating one) and J = 0.0785 for n = 2.5. As illus- 
trated in Fig. 2] the values for a more realistic numerical 
model lie between these curves, as expected. 

Note from Eq. (|2~2"2"|) that s max /^ 2 oc ^~ 2/3 for fi < 
/^max, so that the spin per unit mass squared is un- 
bounded as ii — > 0. 5 Nevertheless, the spin parame- 



ter S'max is bounded, since S„ 



J/J- 



OC /X 



1/3 



4 The mean molecular weight p, is set equal to 2, corresponding 
to helium and heavier elements, which is appropriate for most 
astrophysical white dwarfs. 

5 Eq. |2"3S1 is valid only for u > 0.01 M , but 

Smax//^ continues 

to increase with decreasing [i for equations of state appropriate 
for brown dwarfs and planets. 
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0.6 0.8 

fi(M e ) 

FIG. 1: The maximum spin angular momentum s max 
vs. mass [i for a rigidly rotating white dwarf. We plot curves 
for n = 1.5 and n = 2.5 polytropic approximations us- 
ing Eq. I)2.22|l . together with four points derived using a more 
realistic numerical white dwarf model (Geroyannis and Papa- 
sotiriou [2^|). 



10 



=t 6 




0.2 0.4 



0.6 0.8 

H(Me) 



1.2 1.4 



FIG. 2: s m ax/V vs. mass for a white dwarf. As in Fig. ^ we 
plot curves for n = 1.5 and n — 2.5 polytropes and the numer- 
ical model from [23] . The corresponding spin parameter S max 
is obtained simply by dividing s max //i by the mass M of the 
central black hole. 



the low mass limit. We plot s max //i vs. fi in Fig. |2 
which shows that the maximum value of s max //i is ap- 
proximately 9 Mq (corresponding to a zi = 0.5 Mq white 
dwarf). For a central black hole of mass M = 10 6 M Q , 
we then have 



7 
6.5 

6 
5.5 

5 
4.5 



tides unimportant 



tidal disruption 



0.2 0.4 



0.6 0.8 

A* (Mq) 



1.2 1.4 



FIG. 3: The minimum black hole mass M required not to 
disrupt an inspiraling corotating white dwarf before the last 
stable (prograde) circular orbit around a maximally rotating 
Kerr black hole, as a function of white dwarf mass ti. 



central body overcomes its self-gravity, i.e. 

GM G[i 



(2.24) 



For the white dwarf to be undisrupted down to the hori- 
zon at r = M, we must have M < i? 3 / 2 pT 1 ! 2 , so that 
[using Eq. 1)2.181) ] the minimum mass not to disrupt is 
Afmin oc fi . We could evaluate the proportionality con- 
stant using Eq. (|2.18[) , but we can obtain a more accurate 
result by adopting a constant based on a more realistic 
tidal disruption model. Tables 1 and 2 of |24[ give the 

value of the variable f = — ( — ] , which is approx- 

R\MJ 

imately 2.0 for the white dwarfs of interest here. This 
gives 



M min = 2.0~ 3/2 i? 3 /V 1/2 



(2.25) 



as illustrated in Fig. [3] For a 1.0 Mq white dwarf, 
which (based on [23]) has s max = 8.57 Mq, the cen- 
tral black hole must satisfy M m i n = 8.2 x 10 4 M Q , 
so that the spin parameter S can be no bigger than 
•Smax = s max /(^M min ) = 1.0 x 10~ 4 in order to avoid 
tidal disruption. 



S<S max = ^ = 9x 10- 6 , 



which is small compared to unity. 



(2.23) 



2. Tidal disruption 

We can obtain a higher value of S if the central black 
hole mass is smaller, but it is important to bear in mind 
that such lower-mass black holes may tidally disrupt the 
white dwarf companion, thereby violating a necessary 
condition for the validity of the Papapetrou equations. 
In order of magnitude, a white dwarf orbiting at radius r 
will be disrupted when the tidal acceleration due to the 



3. The S = 1 limit 

We have shown that all physically realistic cases sat- 
isfy S < 1, but we nevertheless consider the limit 
of S — 1 (corresponding to \l = M) in order to inves- 
tigate more thoroughly the dynamics of the Papapetrou 
equations, and to compare our results with 3], which in- 
vestigates the 5=1 limit in detail. The S = 1 limit 
introduces no singularities into the equations of motion, 
and the resulting orbits are valid solutions of the equa- 
tions. On the other hand, in this limit the Papapetrou 
equations are not physically realistic, since they are de- 
rived in the limit of spinning test particles, which must 
satisfy /u -C M. We thus cannot draw reliable results 
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about the behavior of astrophysical systems from the 
5 = 1 limit. 



C. Symmetries and the parameterization of initial 
conditions 

In the approximation represented by the Papapetrou 
equations there is still a constant of the motion associated 
with each Killing vector of the spacetime : 

Q = - it^S^. (2.26) 

[For brevity, we write the constant in terms of the spin 
tensor [Eq. (|2.10() ].] Since Kerr spacetime is sta- 
tionary and axially symmetric, it has the Killing vectors 
£ f = d/dt and £^ = d/d(f>, so the energy E and z angular 
momentum J z are conserved: 

E = -p t + \g t ^vS^ (2.27) 

and 

J*=P*-h94>wS' w . (2.28) 

(We write J z in place of the orbital angular momen- 
tum L z since the spin also contributes to the angular 
momentum of the system.) In contrast to the energy 
and momentum integrals, the Carter constant Q is no 
longer present when the test particle has nonvanishing 
spin |2fij . 

In our problem there are twelve variables, four each 
for position, momentum, and spin. For the purposes of 
finding orbits by numerical integration, we may param- 
eterize the initial conditions by providing r = 0, r, 9, 
4> = 0, p r , E, J z , S, and any two of the spin components. 
The normalization conditions p^p^ = — 1 and S^S^ = S 2 
allow us to eliminate one component each of momentum 
and spin. The constraint p M 5 M = and the integrals of 
the motion then give three equations in three unknowns: 

= ppSugr (2-29) 
E = -Pt + \gt^uS^ (2.30) 
Jz = P<t>- \a^., v S^ (2.31) 

We must solve these equations for the two remaining 
components of p^ and one remaining component of S^. 
In Schwarzschild spacetime these can be solved explicitly 
due to the greater symmetry 0, but in the Kerr case of 
interest here the problem requires numerical root finding. 

We also use a related parameterization method start- 
ing with the Kerr geodesic orbital parameters: eccentric- 
ity e, inclination angle l, and pericenter r p . We derive 
the corresponding energy, angular momentum, and rele- 
vant momenta, and then proceed as above. This method 
is discussed further in Sec. IVII A 31 



III. LYAPUNOV EXPONENTS 

A. General discussion of Lyapunov exponents 

Our method for calculating Lyapunov exponents is 
well-established in the literature of nonlinear dynami- 
cal systems 0, , but accessible treatments are hard 
to find in the physics literature, so we summarize the 
method here. Our discussion is informal and oriented 
toward practical calculation, based on Ref. [T(|; for a 
more formal, rigorous presentation see Eckmann and Ru- 
elle [H. 

First we give an overview of the methods for calculat- 
ing Lyapunov exponents most commonly used in physics. 
Given an initial condition, a set of differential equations 
determines a solution (the flow), which is a curve in 
the phase space. The Lyapunov exponents of the flow 
measure the rate at which nearby trajectories separate. 
As discussed in the introduction, an orbit is chaotic if 
a nearby phase-space trajectory separated by an initial 
distance eo separates exponentially: e(r) = eoe Xr , where 
A is the Lyapunov exponent. 

Implicit in the definition of chaos above is a notion of 
a distance function on the phase space (or, more prop- 
erly, the tangent space to the phase space, as in Eq. i|3.3|) 
below). It is conventional to use a Euclidean metric to 
define such lengths 0,^2' but any positive-definite non- 
degenerate metric will do [2(|. While the magnitude of 
the resulting exponent obviously depends on the particu- 
lar metric used, the signs of the Lyapunov exponents are 
a property of the dynamical system and do not rely on 
any underlying metric structure. We discuss these issues 
further in Sec. HVaI and Sec. I VII PI 

This informal definition of Lyapunov exponents leads 
to a practical method for calculating A: given an initial 
condition, consider a nearby initial condition a distance 
eo away, where eo is "small" , typically 10 -5 -10 -7 of the 
relevant physical scales. (Values of eo much smaller than 
this can result in a loss of numerical precision.) Keep- 
ing track of the deviation vector between the two points 
yields a numerical approximation of A. (It is important 
to rescale the deviation vector if it grows too large, since 
for any bounded phase space flow even a tiny deviation 
can grow to at most the size of the bounded region.) We 
call this approach the deviation vector method. 

There are two primary limitations to the approach out- 
lined above. First, the method yields only the largest 
Lyapunov exponent, which is sufficient to establish the 
presence of chaos but paints a limited picture of the dy- 
namics. Second, the deviation vector approach is most 
appropriate when an analytical expression for the Jaco- 
bian matrix is unknown; by choosing eo small enough 
(and by keeping e(r) small by rescaling if necessary), the 
method essentially takes a numerical derivative. Among 
other complications, the value of the exponent depends 
both on the maximum allowable size e max (the size at 
which the deviation is rescaled) and the initial value eo 
(the size of the deviation after each rescaling). 
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The principal virtue of the deviation vector approach 
compared to the more complicated Jacobian method (dis- 
cussed below) is speed, since it requires solving only the 
equations of motion. (As we discuss in Sec. IIIIB II the 
Jacobian method involves the time-consuming evolution 
of the Jacobian matrix in parallel with the equations of 
motion.) It also provides a valuable way to verify the 
validity of the Jacobian method. 

The Jacobian method is a more thorough and rigor- 
ous approach to the calculation of Lyapunov exponents, 
which makes precise the notion of "infinitesimally" sep- 
arated vectors. The general method proceeds as follows: 
consider a phase space with variables y = {z^} and an 
autonomous set of differential equations 



oIt 



= f(y). 



(3.1) 



(Here we use r instead of t in anticipation of the appli- 
cation of these results to general relativity, where we will 
be using proper time as our time parameter.) If Sy repre- 
sents a small deviation vector, then the distance between 
the two trajectories is 



d(Sy) 
air 



= f (y + Sy) - f (y) = Df • Sy + 0(||<5y|| 2 ), (3.2) 



where Df is the Jacobian matrix [(Df)y = dfi/dx 3 ]. 

We can clarify the notation and make the system easier 
to visualize if we introduce £ as an element of the tangent 
space at y, so that 



dr 



Df 



(3.3) 



which is equivalent to taking the limit ||<5y|| — > 0. We vi- 
sualize £ as a perfectly finite vector (as opposed to an "in- 
finitesimal" ) . Since it lives in the tangent space, not the 
physical phase space, £ can grow arbitrarily large with 
time. This means that instead of the frequent rescal- 
ing required in the deviation vector approach, £ must be 
rescaled only when it grows so large that it approaches 
the floating point limit of the computer. This is a rare oc- 
currence, and in practice the tangent vector almost never 
needs rescaling. 

Although following the evolution of an arbitrary ini- 
tial tangent vector £ yields the largest Lyapunov expo- 
nent, we can do even better by following the evolution 
of a family of n tangent vectors, which allows us to de- 
termine all n Lyapunov exponents. The essence of the 
method is as follows: for a system of differential equa- 
tions with n variables, we consider a set of n vectors that 
lie on a ball in the tangent space. We represent this ball 
using a matrix whose columns are n normalized, linearly 
independent tangent vectors, conventionally taken to be 
orthogonal. This set of orthonormal vectors then spans a 
unit ball in the tangent space. The action of the Jacobian 
matrix, which is a linear operator on the tangent space, is 
to map the ball to an ellipsoid under the time-evolution 
of the flow, as shown in Fig. 0] 




FIG. 4: The Jacobian matrix maps a tangent space ball to 
an ellipsoid. 



For a dynamical system with n degrees of freedom, 
there are n Lyapunov numbers that measure the average 
growth of the n principal axes {rj(r)}™ =1 of the ellipsoid. 
More formally, the Lyapunov numbers Li are given by 



lim h(r)] 1 ^, 



(3.4) 



where Vi (r) is the length of the zth principal axis of the 
ellipsoid. The corresponding Lyapunov exponents are the 
natural logarithms of the Lyapunov numbers, so that 



A, 



lim 

T — >00 



log [r t (r) 



(3.5) 



These limits exist for a broad class of dynamical sys- 
tems 0. 

The principal axes of the tangent space ellipsoid in- 
dicate the directions along which nearby initial condi- 
tions separate or converge, which we may call the Lya- 
punov directions. In particular, consider a principal axis 
that is stretched under the time evolution. Such a vec- 
tor has one component for each dimension (position or 
momentum) in the phase space; a nonzero component 
in any direction indicates an exponential divergence in 
the corresponding coordinate. For example, if a system 
has two spatial coordinates (r, <fi) and corresponding mo- 
menta (prjP</>), then a typical tangent vector will have 
components £ = (^r,^</>,^p r ,^p^)- If the only tangent 
vector with nonzero Lyapunov exponent is, for example, 
£ = -4=(1, 1, 1,0), then nearby initial conditions separate 
exponentially in r, 0, and p r , but nearby values of p^ do 
not separate exponentially. This is potentially relevant 
to the present study since, in the limit of a point test par- 
ticle, the gravitational radiation depends on the spatial 
variables but not the spin. If the principal axes along 
expanding directions have nonzero components only in 
the spin directions, the system could be formally chaotic 
without affecting the gravitational waves. 

In summary, the method for visualizing the Lyapunov 
exponents of a dynamical system is to picture a ball of 
initial conditions — an infinitesimal ball if visualized in 
the phase space, or a unit ball if visualized in the tan- 
gent space — and watch it evolve into an ellipsoid un- 
der the action of the Jacobian matrix. After a suffi- 
ciently long time, the ellipsoid will be greatly deformed, 
stretched out along the expanding directions and com- 
pressed along the contracting directions. The directions 
of the principal axes are the Lyapunov directions, and 
their lengths give the Lyapunov numbers through the re- 
lation L, w [r l (r)] 1 / T . 
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B. Numerical calculation of Lyapunov exponents 

In order to implement a numerical algorithm based on 
the considerations above, we must bear two things in 
mind. First, since the vectors spanning the initial unit 
ball are arbitrary, they will all be stretched in the di- 
rection of the largest exponent: in general every initial 
vector has some nonzero component along the direction 
of greatest stretching, which dominates as r — > oo. In 
order to find the other principal axes, we must periodi- 
cally produce a new orthogonal basis. We will show that 
the Gram-Schmidt procedure is appropriate. Second, the 
lengths of the vectors could potentially overflow or un- 
derflow the machine precision, so we should periodically 
normalize the ellipsoid axes. 



1. The algorithm in detail 

To simplify the notation, we denote the (time- 
dependent) Jacobian matrix Df by J T and the ellipsoid 
(whose columns are the tangent vectors) by U. The al- 
gorithm then proceeds as follows: 

1. Construct a set of n orthonormal vectors (which 
span an n-dimensional ball in the tangent space of 
the flow) . Represent this ball by a matrix U whose 
columns are the tangent vectors £j. 

2. Eq. (|3.3|1 . applied to each tangent vector, implies 
that U satisfies the matrix equation 

^=JrU, (3.6) 

which constitutes a set of linear differential equa- 
tions for the tangent vectors. Since J T depends on 
the values of y, these equations are coupled to our 
system of nonlinear differential equations y = f (y) , 
so they must be solved in parallel with Eq. I|3.1|l . 

3. Choose some time T big enough to allow the ex- 
panding directions to grow but small enough so 
that they are not too big. Numerically integrate 
Eqs. 1|3.1[) and l|3.6l) . and every time T apply the 
Gram-Schmidt orthogonalization procedure. The 
vectors resulting from the Gram-Schmidt proce- 
dure approximate the semiaxes of the evolving el- 
lipsoid. Record the log of the length log [r»(r n )] 
of each vector after each time T, where r„ = nT . 
Finally, normalize the ellipsoid back to a unit ball. 

4. At each time r, the sum 

A^if l g[n(r n )]^i2sMl)l (3 .7) 
r » — ' t 

n—l 

is a numerical estimate for the ith Lyapunov expo- 
nent. 



2. Gram-Schmidt and Lyapunov exponents 

The use of the Gram-Schmidt procedure is crucial to 
extracting all n Lyapunov exponents. Let us briefly 
review this important construction. Given n linearly- 
independent vectors {u^}, the Gram-Schmidt procedure 
constructs n orthogonal vectors {v^} that span the same 
space, given by 

v, = u,-g^ Vi . ( 3. 8) 

To construct the ith orthogonal vector, we take the ith 
vector from the original set and subtract off its projec- 
tions onto the previous i — 1 vectors produced by the 
procedure. 

The use of Gram-Schmidt in dynamics comes from ob- 
serving that the resulting vectors approximate the semi- 
axes of the tangent space ellipsoid. After the first time T, 
all of the vectors point mostly along the principal expand- 
ing direction. We may therefore pick any one as the first 
vector in the Gram-Schmidt algorithm, so choose £i = Ui 
without loss of generality. If we let e, denote unit vec- 
tors along the principal axes and let be the lengths of 
those axes, the dynamics of the system guarantees that 
the first vector Ui satisfies 

ui = nei + r 2 e 2 H « riei = vi 

since ei is the direction of fastest stretching. The second 
vector V2 given by Gram-Schmidt is then 

Ui • Vi 

v 2 = Ui - - — rpr- vi « ui - nei = r 2 e 2 , 

Kir 

with an error of order r^/vx. The procedure proceeds 
iteratively, with each successive Gram-Schmidt step (ap- 
proximately) subtracting off the contribution due to the 
previous semiaxis direction. 

It is important to choose time T long enough to keep 
errors of the form r%jr\ small but short enough to prevent 
numerical under- or overflow. In practice, the method is 
quite robust, and it is easy to find valid choices for the 
time T, as discussed in Sec. IVIII 

IV. RELATIVITY AND PAPAPETROU 
SUBTLETIES 

The algorithm described above is of a general nature, 
designed with a generic dynamical system in mind. The 
Papapetrou equations and the framework of general rel- 
ativity present additional complications. Here we dis- 
cuss some refinements to the algorithm necessary for the 
present case. 

A. Phase space norm 

In the context of general relativistic dynamical sys- 
tems, the meaning of trajectory separation in phase space 
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is somewhat obscured by the time variable. We can skirt 
the issue of trajectories "diverging in time" by using a 
3 + 1 splitting of spacetime, and consider trajectory sep- 
aration in a spacelike hypersurface [27| . This prescription 
reduces properly to the traditional method for classical 
dynamical systems in the nonrelativistic limit. 

In Kerr spacetime, we use the zero angular-momentum 
observers (ZAMOs), and project 4-dimcnsional quanti- 
ties into the ZAMO hypersurface using the projection 
tensor P% = 5^ + U^U V , where U» is the ZAMO 
4-velocity. In this formulation, spatial variables obey 
2+ — * x l = P l £ M and momenta obey — » pi = P^ p^ 
(and similarly for S^) j^]. The relevant norm is then a 
Euclidean distance in the 3-dimcnsional hypersurface. 

We should note that we use the projected norm for 
conceptual clarity, not necessity. The naive use of a Eu- 
clidean norm using unprojected components yields the 
same sign for the exponents, as noted in Sec. IIII Al 
The magnitudes of the resulting exponents are also sim- 
ilar (Sec. ETTDj. 



B. Constraint complications 

Although the Lyapunov algorithm is fairly straight- 
forward to implement for a general dynamical system, 
the constrained nature of the Papapetrou equations adds 
a considerable amount of complexity. The fundamental 
problem is that the tangent vector £ cannot have arbi- 
trary initial components for the Papapetrou system, as it 
can for an unconstrained dynamical system. Each £ must 
correspond to some deviation Sy which is not arbitrary: 
the deviated point y + <5y must satisfy the constraints. 



1. Constraint- satisfying deviations 

Recall that the dynamical variables in the Papapetrou 
equations must satisfy normalization and orthogonality 
constraints (Sec. Ill A|) : p v p v = — 1 (normalized units), 
S v S v — S 2 , andp" 'S„ = 0. To make the constraint condi- 
tion on Sy clearer, let C(y) = represent the constraints 
rearranged so that the right hand side is zero. For exam- 
ple, with y = (t,r, Li 1 (f) 1 p t ,p r ,p f _ l ,p l p, S t , S r , S^, S^,), 6 we 
can write 

C 1 (y)=p"jv + 1, (4.1) 

so that C\ (y) = for a constraint-satisfying y. The other 
constraints are then 

C 2 (y) = S v S„-S 2 (4.2) 



and 

C 3 (y)=p v S v . (4.3) 

A deviation Sy is constraint- satisfying if C(y + Sy) = 
when C(y) = 0. 

We may construct a constraint-satisfying deviation <5y 
as follows. Begin with a 12-dimensional vector y that 
satisfies the constraints. Add a random small deviation 
to eight of its components to form a new vector y'. (We 
need not add a deviation to t; see Sec. IIV B 21 below.) 
Determine the remaining three components of y' using 
the constraints, using the same technique used to set the 
initial conditions. Finally, set Sy = y' — y. The corre- 
sponding £ is then simply <5y/||5y||. 

The prescription above glosses over an important de- 
tail: the inference of tangent vector components from 
the constraints is not unique. Solving the constraint 
equations involves taking square roots in several places, 
so there are a number of sign ambiguities representing 
different solution branches. The implementation of the 
component-inference algorithm must compare each com- 
ponent of y with the corresponding component of y' 
to ensure that they represent solutions from the same 
branches. Enforcing the constraints in this manner, and 
thereby inferring the full tangent vector is especially 
important for the algorithm described in the next section. 



2. A modified Gram-Schmidt algorithm 

A spinning test particle has an apparent twelve de- 
grees of freedom — four each for position, momentum, and 
spin — so a priori there is the potential for twelve nonzero 
exponents. Since the Papapetrou equations have no ex- 
plicit time-dependence, we can eliminate the time degree 
of freedom. The three constraints (momentum and spin 
normalization, and momentum-spin orthogonality) fur- 
ther reduce the number of degrees of freedom by three. 
We are left finally with eight degrees of freedom. 

Eliminating the four spurious degrees of freedom from 
the tangent vectors presents a formidable obstacle to 
the implementation of the phase space ellipsoid method 
described in Sec. IIII B II The crux of the dilemma 
is that the axes of the ellipsoid must be orthogonal, 
but must also correspond to constraint-satisfying devia- 
tion vectors — mutually exclusive conditions. Solving this 
problem requires a modification of the Gram-Schmidt al- 
gorithm: 



1. Instead of a 12 x 12 ball (i.e., n = 12 in the original 
algorithm), consider an 8 x 8 ball by choosing to 
eliminate the t, p t , p^, and St components. The 
time component £t of each tangent vector is irrele- 
vant since nothing in the problem is explicitly time 
dependent; the first column of the Jacobian ma- 
trix is zero, so £t is not necessary to determine the 



Recall that we write the equations of motion in terms of fj, = 
cos 6. 
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time-evolution. 7 The other three components are 
determined by the constraints as described above. 

2. Given eight initial random tangent vectors, apply 
the Gram-Schmidt process to form an 8 x 8 ball. 
For each vector, determine the three missing com- 
ponents using the constraints, and then evolve the 
system using 



as before. (Now U represents a 12 x 8 matrix in- 
stead of a 12 x 12 ball.) 

3. At each time T, extract the relevant eight compo- 
nents from each vector to form a new 8x8 ellipsoid, 
apply Gram-Schmidt, and then fill in the missing 
components using the constraints, yielding again a 
12x8 matrix. The projected norms of the eight tan- 
gent vectors contribute to the running sums for the 
Lyapunov exponents as in the original algorithm. 

The algorithm above yields eight Lyapunov exponents for 
the Papapetrou system of equations. 

In order to implement this algorithm, we must have a 
method for constructing a full tangent vector £ from an 
eight-component vector £. The method is as follows: 

1. Let y' = y + e£ for a suitable choice of e. 

2. Fill in the missing components of y' using the con- 
straints to form y', taking care that y and y' have 
the same constraint branches. 

y' — y 

3. Infer the full tangent vector using £ = . 

e 

This technique depends on the choice of e, and fails 
when e is too small or too large. Using the techniques 
discussed in the next section to calibrate the system, we 
find that e « 10 -5 -10 -6 works well in practice. 



the constraint subtleties. The complexity of the main Ja- 
cobian approach involves the competing requirements of 
Gram-Schmidt orthogonality and constraint satisfaction, 
but in the case of only one vector these difficulties vanish. 
Since the equations of motion preserve the constraints, 
an initial constraint-satisfying tangent vector retains this 
property throughout the integration. Thus, we begin 
with a vector constructed as in Sec. lIVBTl and evolve it 
(without rescaling) along with the equations of motion. 
Other than the requirement of constraint satisfaction, its 
initial components are arbitrary, so it evolves in the di- 
rection of largest stretching and eventually points in the 
largest Lyapunov direction. The logarithm of its pro- 
jected norm then contributes to the sum for the largest 
Lyapunov exponent. 

Second, we can implement a deviation vector ap- 
proach as described in Sec. IIII Al Given an initial con- 
dition yo, we construct a nearby initial condition yQ as 
in Sec. II V B II and then evolve them both forward. In 
principle, an approximation for the largest Lyapunov ex- 

, . i , ( \W-y\\ \ i , ( WW \ t 

ponent is then - log — = - log j- — - . In 

t Vllyo-yolly r \Wo\\J 

practice (for chaotic systems) the method saturates: for 
a given initial deviation, say ||<5yo|| ~ 10~ 6 , once the 
initial conditions have diverged by a factor of ~10 6 the 
method breaks down. 8 (The traditional solution to the 
saturation problem is to rescale the deviation before it 
saturates, but such a rescaling in this case violates the 
constraints.) Despite its limitations, this unrescaled de- 
viation vector technique is valuable, since it tracks the 
correct solution until the saturation limit is reached, and 
avoids the subtleties associated with the constraints. 

With these two techniques in hand, we have a powerful 
method for verifying that the largest Lyapunov exponent 
produced by the Gram-Schmidt method is correct. This, 
in turn, gives us confidence that the other Lyapunov ex- 
ponents produced by the main algorithm are meaningful 
as well. 



3. Two rigorous techniques 

It should be clear from the discussion above that ex- 
tracting all eight Lyapunov exponents is difficult, and 
in practice the techniques are finicky, depending (among 
other things) on the choice of e as described in Sec. lIVB"2l 
above. How, then, can we be confident that the re- 
sults make sense? Fortunately, there are two techniques 
that give rigorous Lyapunov exponents by managing to 
sidestep the constraint complexities entirely. 

First, it is always possible to calculate the single largest 
exponent using the Jacobian method without considering 



V. IMPLEMENTATION DETAILS 
A. Some numerical comments 

Finally, we discuss some specialized issues related to in- 
tegrating the Papapetrou equations on a computer. The 
primary subjects are the formulation of the equations, 
optimization techniques, and error checking. 

Our choice to write the Papapetrou equations using 
the spin vector is motivated partially by numerical con- 
siderations. The spin vector approach has nice properties 



7 Also, the time piece is discarded in the projected norm formalism 
in any case (Sec. HV^l . 



8 This underscores the point that chaos is essentially a local phe- 
nomenon. Any unrescaled deviation vector approach must satu- 
rate, since no bounded system can have trajectories that diverge 
for arbitrarily long times. 
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compared to the tensor approach as S — * 0. Comparing 
their covariant derivatives is instructive: 

VaS^ = p»v v = 2p['V / l. 

Though simpler in form, the derivative of S^" has unfor- 
tunate numerical properties for small S, since in the limit 
S — > we have p^ — > the difference — p^w^ goes 
to zero in principle but in practice is plagued by numer- 
ical roundoff errors. Since S >C 1 is the most physically 
interesting limit, the vector approach is more convenient 
for our purposes. 

Calculating the many tensors and derivatives which 
go into the Papapetrou equations and the correspond- 
ing Jacobian matrix is a considerable task. As a first 
step, we use GRTensor for Maple to calculate all rele- 
vant quantities, and we use Maple's optimized C output 
to create C code automatically. Due to the symmetries 
of the Riemann tensor and the metric, many terms are 
identically zero, which significantly reduces the number 
of required operations. For example, in order to calcu- 
late R* a p jS S a v^p 1 Ss we need four loops, which consti- 
tutes 4 4 = 256 evaluations, but in fact R* a ^ lS has only 
80 nonzero components. Performing loop unrolling by 
writing these terms to an optimized derivatives file con- 
sisting of explicit sums speeds up calculation by an order 
of magnitude compared to nested for loops. 

Another optimization involves the choice of coordi- 
nates used in the metric, which has significant conse- 
quences for the size of the tensor files and the num- 
ber of floating point operations required. Simply using 
/i = cos 9 in the Kerr metric reduces the size of the Rie- 
mann derivatives by at least a factor of two. 9 Since these 
derivatives are the bottleneck in the calculation of the Ja- 
cobian matrix, we can get more than a 50% improvement 
in performance with even this simple variable transfor- 
mation. 

All integrations were performed using a Bulirsch-Stoer 
integrator adapted from Numerical Recipes |28j . Occa- 
sional checks with a fifth-order Runge-Kutta integrator 
were in agreement. We verified the Papapetrou integra- 
tion by checking errors in the constraints and conserved 
quantities; for an orbit such as that shown in Fig. all 
errors are at the 10 -11 level after r = 10 5 M. 

As should be clear from Sec. IV Bl below, the Jacobian 
matrix of the Papapetrou equations has a large number of 
terms, and it is essential to verify its correctness by using 
a diagnostic that compares Df • Sy with the difference 



9 Warning: This variable substitution changes the handedness of 
the coordinate system, since the unit vector fi points opposite 
to 8. This in turn introduces an extra minus sign in the Levi- 
Civita tensor e a P~< s , which appears many times in the Papa- 
petrou equations and the corresponding conserved quantities. 
The author discovered this subtlety the hard way. 



f(y + <fy) — f (y) for a suitable constraint-satisfying Sy. 
It is not sufficient for the difference merely to be small: 
we must calculate the quantity f (y + Sy) — f (y) — Df • Sy 
for several values of Sy and verify that each component 
scales as ||<5y|| 2 - An early implementation of the Jacobian 
matrix, which gave nearly identical results for f {y+Sy) — 
f (y) and Df • Sy, nevertheless had an undetected 0(S 2 ) 
error. The unrescaled deviation vector approach showed 
a discrepancy with the Jacobian method, 10 which showed 
spurious chaotic behavior. The ||<5y|| 2 scaling method 
described above eventually diagnosed the problem, which 
resulted from a missing term in dS^/dSv (Sec IV B} . 

B. The Jacobian matrix 

For reference, we write out explicit equations for part 
of the Jacobian matrix of the Papapetrou equations. 

The Jacobian matrix of a system of differential equa- 
tions, specialized to the case at hand, is as follows: 

/ dx^ dx^ dx* 1 \ 

~dx» ~dp~ v dS^ 

dp^ dp^ dp^ 

dx v dp v dS v 1 ' 

dS^ dS^ dS^ 
\ dx v dp v dS v J 

Once we calculate dx^/dx 1 " = all the other deriva- 
tives can be expressed in terms of the derivatives of v M , 
the tensors and connection coefficients, and Kronecker 
5s. 

Written out in full, and the Papapetrou equations are 
as follows: 

±" = «" (5.2) 

- -R% a0 vy a Sp + T a ^p a vP (5.3) 

- - Pl ,(R*y s S a v%S s ) +T%S a vP (5.4) 

We measure r and r in units of M (the mass of the central 
body), p^ in units of the particle rest mass /z, and in 
terms of the product fiM . The over dot is an ordinary 
derivative with respect to proper time: x = dx/dr. 

The unusual placement of indices on R* is motivated 
by the form of the Jacobian matrix. The index placement 
shown above brings the equations into a form where the 
indices on p M and are always lowered, which simplifies 
the Jacobian matrix since (for example) dp^/dx^ = 0. 
Otherwise the Jacobian matrix is unnecessarily compli- 
cated; for example, if p^ appeared anywhere on the right 
hand side then we would have dp^ jdx v ^ 0, which would 
contribute to J T . 



This illustrates the value of calculating the Lyapunov exponents 
using two different methods. 
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As discussed in Sec III Al the supplementary condition 
p^S^ = [Eq. (|21J] leads to the equation for in 
terms of p^ : 



with 



where 



and 



vr 



, R*^S a p p S 1 



(5.5) 
(5.6) 
(5.7) 



N is a normalization constant fixed by v^v^ = — 1. 

The calculation of the partial derivatives i M in 
Eq. 1)5.1(1 proceeds as follows. From the relation for 
= jVw^, we have 



dx v 



— „M — AJ^f 1 



Now, v% = p% + w» tU = Pa g a ^ - *R*^ u S a ppS 1 , 
so the first term is easy. The second term is trickier: 
from the expression for v M , we have that — 1 = v^v^ — 
N 2 ij^p^ + 2w' J 'p ft + #w M ) = N 2 (-l + 2w^p^ + w^Wfj), 
so we have 

N = (1 - 2w% - w^w^)- 1/2 . 
Differentiating gives 

N, v = N 5 (p a w a u + w%w a + \w a wPg a ^ v ) 

where we have re-labeled the dummy index (/i — > a). 
Summing the various terms, we have 

<„ = N[p a g a ^ + w^ 

+ v^{v a w% + \NvPwPg af ,, v )] . (5.8) 

The expression for dx^ jdp v is similar to but it is 
simpler because the derivative of the metric with respect 
to the momentum is zero. As before, we use the product 
rule: 

g v » Qv" ON _ u 

Tj = "t ' T ' 

op v ap v op v 
The first term requires 

W = l^ + dwP_ = _ * R * m »0 SaS = g »v + 
Op,, Op v Op v 

Note that W^ u is symmetric. The second term requires 
|^ = N a {W au p a +w a S a u +W av w a ) = N\w»+v a W a »). 
Summing the terms gives 



dpu 



W» v = -*R*^S a Sf3. (5.10) 



Finally, we calculate dx^/dS^. With 

™- = -S a p p {*R*^ af3y - *R*^ a P) = v^, 



and 



we have 



ON 

— &< = Nv»v a V au , 



dp v 



Nig" 1 ' + + Nv^w") + Nv^v a W av (5.9) 



— = NV^ + Nv"v a V au . (5.11) 

We calculate the derivatives of p M and using u%, 
the product rule, and the derivatives of the various ten- 
sors in the problem. The full results appear in the ap- 
pendix. 



VI. INTEGRABILITY AND CHAOS 

A. Phase space and constants of the motion 

Having laid the foundation for the numerical calcula- 
tion of Lyapunov exponents, we now discuss some general 
aspects of dynamical systems relevant to our study. A dy- 
namical system with n coordinates has a 2n dimensional 
phase space, typically consisting of generalized positions 
and their corresponding conjugate momenta. Motion in 
the phase space is arbitrary in general, but when there 
are integrals of the motion then the flow is confined to 
a surface on which the integral is constant. This can be 
seen most easily by transforming to angle-action coor- 
dinates, where the surface is an invariant (multidimen- 
sional) torus. 

A system with n coordinates and n constants of the 
motion is integrable and cannot have chaos (though the 
motion can still be quasiperiodic or exhibit other compli- 
cated behavior). For example, we can consider geodesic 
orbits around a Kerr black hole to have eight degrees 
of freedom (n = 4) and four constants of the motion — 
particle rest mass /x, energy E, axial or z angular mo- 
mentum L z , and Carter constant Q — which are enough 
to integrate the equations of motion explicitly. Alter- 
natively, we may look at Kerr spacetime as having a 6- 
dimensional phase space by eliminating time (which is 
simply a re-parameterization of the proper time) and us- 
ing rest mass conservation to eliminate one momentum 
coordinate. Then the three integrals E, L z , and Q are 
sufficient to integrate the motion. (In practice, we allow 
all four momenta to evolve freely; the normalization is 
then a constraint which can be checked for consistency 
at the end of the integration.) 
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In the case of a spinning test particle, the extra spin 
degrees of freedom create the possibility for chaotic be- 
havior. Moreover, since Q is not conserved in the case of 
nonzero spin, even without the extra spin degrees of free- 
dom the potential for chaos would exist. Kerr spacetime 
has just enough constants to make the system integrable; 
losing Q reduces the number of analytic integrals below 
the critical level required to guarantee integr ability. 11 



B. Hamiltonian systems 

1, Lyapunov exponents for Hamiltonian flows 

The phase space flow of Hamiltonian systems is con- 
strained by more than the integrals of the motion. In 
particular, the Lyapunov exponents of a Hamiltonian sys- 
tem come in pairs ±A; i.e., if A is a Lyapunov exponent 
then so is —A j2^|. Geometrically, this means that if one 
semimajor axis of the phase-space ellipsoid stretches an 
amount e Ar = L, another axis must shrink by an amount 
e~ Ar = 1/L. One consequence of this property is that 
the product of the lengths of the axes is 1. Since the el- 
lipsoid volume is proportional to this invariant product, 
Liouville's theorem on the conservation of phase space 
volume follows as a corollary. 

The ±A property of Hamiltonian flows results from the 
symplectic nature of the Jacobian matrix for Hamilto- 
nian dynamical systems. 12 But a naive analysis of the 
Jacobian matrix of the Papapetrou equations shows that 
it is not symplectic in the canonical sense. Neverthe- 
less, the Papapetrou equations can be derived from a 
Lagrangian 30] , and can be cast in Hamiltonian form by 
use of a free Hamiltonian with added constraints (follow- 
ing the method of Dirac j^H as discussed in |32|). As 
a consequence, we could in principle find coordinates in 
which the Jacobian matrix is symplectic with respect to 
the canonical symplectic matrix. Fortunately, this is an 
unnecessary complication, since the underlying dynamics 
are independent of the coordinates. 



2. Exponents for spinning test particles 

As discussed in Sec. UVB^l the lack of explicit time de- 
pendence independence and the three constraints reduce 
the degrees of freedom from twelve to eight, which leaves 
the possibility of eight nonzero Lyapunov exponents. The 



11 It is possible that deformations of Kerr geometry that destroy Q 
nevertheless possess a numerical integral that preserves integra- 
bility, in analogy with some galactic potentials 121-1 . but the loss 
of Q certainly ends the guarantee of integrability. 

12 A matrix S is symplectic with respect to the canonical symplectic 

matrix J if S T JS = J, where J = ( ^ j an< ^ ^ ' S n X n 
identity matrix. 



phase space flow is further constrained by the constants 
of the motion, energy and z angular momentum; corre- 
sponding to each constant should be a zero Lyapunov ex- 
ponent, since trajectories that start on an invariant torus 
must remain there. This leaves six exponents potentially 
nonzero. Since the exponents must come in pairs ±A, 
there should be at most three independent nonzero ex- 
ponents. 



VII. RESULTS 

First we give results for the dynamics of the Pa- 
papetrou equations in the extreme (and unphysical) 
limit S = 1, which represents a violation of the test- 
particle approximation but is still mathematically well- 
defined. We find the presence of chaotic orbits (in agree- 
ment with 0]). We next examine the effects of varying S, 
including the limit S <C 1. Finally, we investigate more 
thoroughly the dynamics for physically realistic spins. 



A. Chaos for S = 1 

1. Maximally spinning Kerr spacetime 

In a background spacetime of a maximally spinning 
Kerr black hole (a = 1) there are unambiguous pos- 
itive Lyapunov exponents for a range of physical pa- 
rameters when 5 = 1. We show a typical orbit that 
produces nonzero Lyapunov exponents in Fig. [BJ The 
orbit has energy E — 0.8837 /i, z angular momentum 
J z = 2.0667 /i-M, and the radius ranges from pericenter 
r p = 1.7 M to apocenter r a — 6.7 M. The Lyapunov 
integrations typically run for 10 4 M, which corresponds 
approximately to 400 0-orbital periods. 

We can illustrate the presence of a chaotic orbit by 
plotting the natural logarithm of the ith ellipsoid axis 
log [rj(r)] vs. t [Eq. 13.7fl ]. so that the slope is the Lya- 
punov exponent, as shown in Fig.[7] 13 There appear to 
be two nonzero Lyapunov exponents; the third largest 
exponent is consistent with zero, as shown in Fig. El The 
reflection symmetry of the figure is a consequence of the 
exponent pairing: for each line with slope A, there is a 
second line with slope —A. 

The main plot in Fig. Ufa) is generated by the modi- 
fied Gram-Schmidt algorithm fSec. IIVB 2(1 . Recall that 
this method depends on the value of e used to infer the 
tangent vector; we find a valid e by calibrating it using 
the rigorous Jacobian method, which must yield an ex- 
ponent that matches the largest exponent from the mod- 



It is traditional to plot log [ri(r)]/r, which converges to the Lya- 
punov exponent as r — ► oo, but it is much easier to identify 
the linear growth of log [^(r)] than to identify the convergence 
of log [ri(r)]/r. The ±A property is also clearer on such plots. 
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FIG. 6: The orbit of a maximally spinning (S = 1) particle in maximal Kerr spacetime, for E = 0.8837 [i and J z = 2.0667 \iM 
(the same values as in Fig. El above). The spin has initial values of 5"" = S M = 0.1, corresponding to an initial angle of 54° 
with respect to the vertical in the particle's rest frame. As in Fig. [13 we plot y vs. x in (a) and z vs. p in (b). The spin causes 
significant deviations from geodesic orbits. 



ified Gram-Schmidt method. The plot in Fig. [7Ja) rep- 
resents the case e = 10~ 6 ; it is apparent that the two 
methods agree closely. The unrescaled deviation vector 
method provides an additional check on the validity of 
the largest exponent, as shown in Fig. Efb). As expected, 
the unrescaled approach closely tracks the full Jacobian 
approach until it saturates. 

The numerical values of the exponents are shown in 
Table [I] The ±A property is best satisfied by ±A ma x, 
the exponents with the largest absolute value. The ex- 
ponents are least-squares fits to the data, with approx- 
imate standard errors of 1%. These errors are not par- 
ticularly meaningful since the exponents themselves can 
vary by ^10% depending on the initial direction of the 
deviation vector. Moreover, even exponents that appear 
nonzero may be indistinguishable from zero in the sense 
of Fig. |H1 for such exponents a error on the fit is 

meaningless. 

For initial conditions considered in Fig. and other 



orbits in the strongly relativistic region near the horizon, 
typical largest Lyapunov exponents are on the order of 
a few x 10~ 3 /M. For the particular case illustrated in 
Fig. El we have A max w 5 x 10~ 3 M _1 , which implies an 
e-folding timescale of T\ = 1/A « 2 X 10 2 A/. This is 
strongly chaotic, with a significant divergence in approx- 
imately eight 0-orbital periods. 

Based on integrations in the case of zero spin, which 
corresponds to no chaos (Lyapunov exponents all zero), 
we can determine how quickly the exponents approach 
zero numerically. 14 Fig. compares the four apparently 
positive exponents with a known zero exponent. Only 
two of the four exponents are unambiguously distinguish- 



As noted in the introduction, it is possible for integrable but 
unstable orbits to have positive Lyapunov exponents. We avoid 
this issue by choosing a baseline orbit that is not unstable. 
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FIG. 7: Natural logarithms of the phase space ellipsoid axes vs. proper time in Kerr spacetime with S = 1. The slopes of 
the lines are the Lyapunov exponents; the largest exponent is approximately A max = 5 x 10 -3 M -1 . The initial conditions 
are the same as in Fig. El an d one point is recorded at each time T = 100 M (Sec. 1111 B It . (a) Full Gram-Schmidt Jacobian 
method (light) with rigorous Jacobian method (dark). The full GS method is rescaled at each time T according the algorithm 
in Sec. 1111 Bl while the rigorous Jacobian method is unrescaled. The two methods agree closely on the value of the largest 
Lyapunov exponent, (b) Rigorous Jacobian method compared to unrescaled deviation vector method. Note that the latter 
method, which started with a deviation of size 10 -7 , saturates at ~16. This corresponds to a growth of e 16 ~ 9 x 10 6 , which 
means that the separation has grown to a size of order unity. 
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FIG. 8: Ellipsoid axis lengths from the upper half of Fig. I7f a) 
(light), compared to an integration with zero spin and hence 
zero Lyapunov exponent (dark). Only two of the four lines 
represent exponents distinguishable from zero. 



TABLE I: Lyapunov exponents in Kerr spacetime in units 
of 10 _3 M _1 , using a least squares fit. The exponents cor- 
respond to the semimajor axis evolution shown in Fig. |7Ja). 
As is typical with the Gram-Schmidt Jacobian method, the 
exponents with the largest magnitudes are determined most 
accurately, and thus show the ±A property most clearly. The 
standard errors on the fit are ~1% for each exponent, but 
these errors are dominated by two systematic errors: (1) the 
variation due to different choices of initial (random) tangent 
vectors; (2) nonzero numerical values even for exponents that 
converge to zero eventually. In particular, the four smallest 
exponents (in absolute value) are indistinguishable from zero 
(see Fig. |H). 



+A 


5.5 


1.5 


0.56 


0.25 


—A 


5.3 


1.6 


0.76 


0.072 



able from zero, consistent with the argument in Sec. I VI 51 
that there should be at most three independent nonzero 
exponents. 

Finally, we note that the components of the direction of 
largest stretching are all nonzero in general. The chaos 
is not confined to the spin variables alone, but rather 
mixes all directions. This indicates that chaos could in 
principle manifest itself in the gravitational waves from 
extreme mass-ratio binaries — but see Sec. IVIlTl below. 



2. Schwarzschild spacetime revisited 

We now reconsider the case of a spin S = 1 particle 
in Schwarzschild spacetime, as investigated in Ref. 
Fig. [3] shows an orbit similar to a chaotic orbit consid- 



ered there (Fig. 4(d) in 3]). A plot of log fa (r)] vs. r 
(Fig. II Op t shows behavior similar to that in Fig. In 
particular, the ±A symmetry is present, apparently with 
two positive exponents. (The other lines are indistin- 
guishable from zero, again using S = orbits as a base- 
line.) The largest exponent of 1.5 x 10 _3 M _1 agrees 
closely with the value from Ref. |3j, which reported an 
exponent of ^2 x 10~ 3 M" 1 for a similar orbit. (This 
agreement is somewhat surprising, since appears not 
to have taken the constrained nature of the deviation vec- 
tors into account. Luckily, the exponents are robust, and 
even unconstrained deviation vectors give nearly correct 
results.) 
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FIG. 9: The orbit of a maximally spinning (S = 1) test particle in Schwarzschild spacetime for E = 0.94738162 fi and 
J z = 4.0 [iM As before, we plot (a) y vs. x and (b) z vs. p — \J x 2 + y 2 . 




5000 10000 15000 20000 5000 10000 15000 20000 

T (M) T (M) 

(a) (b) 



FIG. 10: Natural logarithms of the phase space ellipsoid axes vs. proper time in Schwarzschild spacetime with 5=1. The 
largest exponent is A max ~ 1.2 x 10 -3 M _1 . The initial conditions are the same as in Fig. |U] (a) Full Gram-Schmidt Jacobian 
method (light) with rigorous Jacobian method (dark), (b) Rigorous Jacobian method compared to unrescaled deviation vector 
method. As in Fig.|7Jb), the unrescaled method eventually saturates. 



TABLE II: Lyapunov exponents in Schwarzschild spacetime 
in units of 10 _3 A/ _1 , using a least squares fit. The expo- 
nents correspond to the semimajor axis evolution shown in 
Fig. llUf a'l. which is similar to the orbit in Fig. 4(d) of Ref. Q. 
As with the Kerr case (Table[IJ, the standard errors on the fit 
are ~1% for each exponent, and the same caveats apply. The 
four smallest exponents (in absolute vale) are indistinguish- 
able from zero in the sense of Fig. |H| 



+A 


1.2 


0.67 


0.21 


0.0063 


—A 


1.5 


0.57 


0.10 


0.00023 



3. Kerr and Schwarzschild compared 

The Kerr and Schwarzschild Lyapunov exponents of 
the previous two sections are not all that different; both 
are 10 -2 -10 -3 M _1 in order of magnitude. Nevertheless, 
the two systems prove to be quite different: chaotic orbits 



are easy to find in Kerr spacetime for nearly any initial 
condition that explores the strongly relativistic region 
near the horizon, whereas nearly all analogous orbits in 
Schwarzschild spacetime are not chaotic. 

Fig. ^] compares Kerr and Schwarzschild orbits with 
the same inclination angle i = 10° and eccentricity e = 
0.5 but varying pericenters r p . (Details of this parame- 
terization method, mentioned above in Sec. Ill CI appear 
in [? ].) We insure that the systems are analogous by 
using orbits of S = 1 particles with the same values of 
r p /r ms , where r ms is the radius of the marginally sta- 
ble orbit in the corresponding S — (geodesic) case. 
We use a Kerr geodesic integrator developed by Scott 
Hughes [33| to find r ms , which is the smallest pericenter 
that still yields a stable orbit. For the values of t and e 
considered, r ms = 1.0 M for Kerr and r ms = 4.67 M for 
Schwarzschild. 

It is evident from Fig. ^2 that the Kerr orbits are 
chaotic for a broad range of pericenters, with the maxi- 
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FIG. 11: Comparison of maximally spinning (S = 1) particle 
orbits in Kerr (dark) and Schwarzschild (light). We plot the 
largest Lyapunov exponent versus pericenter (normalized by 
the marginally stable radius). The Kerr initial conditions 
for the innermost orbits are essentially as in Fig. |S| The 
Schwarzschild orbits are identical to their Kerr counterparts 
in inclination (10°) and eccentricity (e = 0.5) but have the 
Kerr parameter a set to zero. The Schwarzschild orbits have 
exponents indistinguishable from zero over the entire range of 
parameters. 



mum Lyapunov A max generally decreasing as the pericen- 
ter increases. In contrast, the Schwarzschild orbits are 
not chaotic anywhere over the entire range of valid ini- 
tial conditions. In fact, we are unable to find any chaotic 
orbits in Schwarzschild spacetime other than the types 
identified by Suzuki and Maeda which were excep- 
tional cases of orbits on the edge of a generalized effec- 
tive potential. In Kerr, on the other hand, chaotic orbits 
appear to be the rule for pericenters near r ms . 



B. Dependence on S 

Since chaos must disappear as S — » 0, we expect to 
see the largest Lyapunov exponent approach zero in this 
limit. This is indeed the case: in Fig. ^] which shows 
the variation of A max with S for two different orbits, we 
see that the chaos unambiguously present when 5 = 1 
is not present for smaller values of S. In particular, the 
largest Lyapunov exponent is indistinguishable from zero 
over the entire range 10 -6 < S < 10 _1 . (The far left of 
the plots have data points for each decade in this range.) 

Although the strength of the chaos generally decreases 
with S, one remarkable feature of Fig. H^T a) is the return 
of chaotic orbits between S ~ 0.25 and 0.4 after their 
disappearance at S ~ 0.5. The effect is qualitatively clear 
in Fig. E| This chaotic "bump" in A max vs. S illustrates 
an important theme in nonlinear dynamical systems: the 
only way to determine whether an orbit is chaotic is to do 
the calculation. Though we certainly expect the strength 
of chaos to be smaller for S < 1 than for S « 1, it is 
impossible, in general, to determine a priori whether a 
particular set of parameters will lead to chaotic behavior. 



C. Physically realistic spins 

The Papapetrou equations are only realistic in the 
test-particle limit, so physically realistic spins must sat- 
isfy S <C 1 (Sec. Ill B^ . This corresponds to likely sources 
of gravitational waves for LISA [3 US H||, e.g., maxi- 
mally spinning fi = 10 Mq black holes spiraling into su- 
permassive M = 10 6 M© Kerr black holes, which have 
spin parameters of S = fi/M = 10~ 5 . Because of their 
likely importance as emitters of gravitational waves, it is 
essential to understand the dynamics of such systems. 

1. Vanishing Lyapunov exponents 

We would like to be able to make a definitive statement 
about the presence or absence of chaos for "small" spins, 
e.g., values of S in the range 10 _2 -10 -6 . Unfortunately, 
when determining Lyapunov exponents numerically, it 
is impossible to conclude definitively that an orbit is or 
is not chaotic, since to do so would require an infinite- 
time integration. On the other hand, for suspected non- 
chaotic orbits, we can provide an approximate bound on 
the e-folding timescale. 

The numerical values of exponents suspected to be zero 
depend strongly on the time of the integration. For exam- 
ple, for values of S in the range 10~ 2 < S < 10~ 6 , the ex- 
ponent in Fig.ll2lappears to be A max « 5x 10~ 4 M^ 1 , but 
this plot represents an integration time of only 10 4 M. 
Longer integration times give correspondingly smaller es- 
timates for the suspected zero exponents (Fig. I14f) . For 
the system shown in Fig. ^] an integration of 10 7 M 
yields an estimate of A max « 3.0 x 10 -7 M _1 for all spins 
in the range 10 -2 < S < 10~ 7 . In this case, the relevant 
Lyapunov timescales are at least 3 x 10 6 M, and are prob- 
ably much longer; the size of the bound is limited only 
by our patience and computer budget. It seems highly 
likely that such orbits are not chaotic. 

2. Spin-induced phase differences 

Even if their Lyapunov exponents are zero, small spins 
affect the relative phase of the orbits, and since phase 
differences accumulate secularly 37], the spin can still 
affect the gravitational wave signal. It is therefore use- 
ful to have a sense of the orders of magnitude of such 
spin-induced phase-shifts. Tables lITTl and lTVI show typical 
values for the phase difference Acf> = <^ gC odesic — '/'spin for 
S = 10~ 5 , where the geodesic and spin systems start with 
the same initial 4- velocity v M . The most useful quantity 
in practice is the phase shift as measured by observers at 
infinity, so we integrate in terms of the Boyer-Lindquist 
coordinate time t in place of r. (This involves multi- 
plying the differential equations by dr/dt at each time 
step.) As is apparent from the tables, the phase shifts 
range broadly, from 10 _1 to 10~ 5 radians after 2000 ra- 
dial orbital periods, but tend to decrease in magnitude 
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FIG. 12: Variation of largest Lyapunov exponent vs. S. (a) The spin S — 1 initial conditions are the same as in Fig. |S| 
(b) Another S = 1 case with a different inclination angle (20°) and pericenter (2.5 M). As the spin decreases, we hold fixed the 
Kerr orbital parameters: inclination angle, eccentricity, and pericenter. Note that in (a) the chaos disappears below S ~ 0.5, but 
returns in a region centered on S ~ 0.3. The horizontal line in both plots is the value of A max calculated for the baseline S = 
orbit. In both (a) and (b) the Lyapunov exponent is indistinguishable from zero for physically realistic spins. 



TABLE III: Phase shifts A<f) — gco dcsic — 0s P m in radians as a function of orbital inclination angle i and pericenter r v for 
a — 0.5 and S = 10~ J . Inclination angle l = 0° is prograde equatorial and t = 180° is retrograde equatorial. The geodesic 
orbits and their corresponding spin orbits start with the same initial 4- velocity u M , and the integrations are performed using 
Boyer-Lindquist coordinate time t, with i max ~ (2000 times the average radial orbital period). The pericenters are scaled 
by the marginally stable radius r ms , and we start at r p /r ms = 1.5 to guarantee the existence of valid initial conditions for 
the non-geodesic orbit. The spin has fixed initial values of 5* r = 5* M = 0.1 S (with hats indicating an orthonormal basis), 
corresponding to initial angles of 9° to 30° with respect to the vertical in the particle's rest frame, increasing with decreasing 
pericenter. 
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FIG. 13: Two orbits from the "bump" in Fig. H^Tal. The 
S — 0.4 orbit (light) is not chaotic, but the S = 0.3 orbit 
(dark) is chaotic, despite having a smaller value of the spin. 



with increasing inclination angle or pericenter. 

Ref. 36] shows that the number of orbital periods in 
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FIG. 14: The variation of the dimensionless quantity A max M 
with final integration time Tf for spin parameter S in the 
range 10 -2 < S < 10 -6 . From top to bottom, the total 
integration time varies from 10 4 M to 10 7 M. It is likely that 
the true Lyapunov exponent is zero. 
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TABLE IV: Phase shifts A(f> = </> gC odcsic — <f>apin in radians as a function of orbital inclination angle i and pericenter r v for 
a — 1 and S = 10~ 5 . As in Table ITTT1 the pericenters are scaled by r mB , and the spin has fixed initial values of S r — = 0.1 S 
(corresponding in this case to initial angles of 28° to 61°, again decreasing with increasing pericenter). 
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FIG. 15: The natural logarithm of the largest ellipsoid axis 
vs. time for proper time r (dark) and coordinate time t 
(light). The exponents are clearly different, but the Lyapunov 
timescales t> = 1/A T and t\ = 1/A t are related by Eq. 17.11 . 



a full inspiral from r « 4M to the final plunge is N ~ 
which is 10 5 for the systems in Tables ITTT1 and llVl 
Since the table represents values of A<f> for 2000 times 
the average radial orbital period, this means that the 
total phase shift during the inspiral is 50 Arabic- For 
a 10° inclination angle the total phase shift is on the 
order of a tenth of a radian to a radian. Slightly more 
realistic values of the number of orbits can be obtained 
using Fig. 2 in |36l |. which gives N ~ 2 x 10 4 orbital 
periods from r — AM to the plunge at r r* M for a = 
0.998, i = 10°, and M/(i = 10 5 . Since the orbit spends 
most of its time between 4r ms and 2r ms , interpolating 
in Table Hvl gives A0 to tai ~ 10 X A0 r=3 . o = 2 x 10" 2 . 
This is only a rough estimate, since the orbits in |3fij | are 
circular, while the orbits we consider are eccentric. 



D. Comments on time, rescaling, and norms 

In this paper, we have elected to use r as the time 
parameter, a rescaling time T of 100 M, and a projected 
norm (Scc. lIV A"|) . Here we discuss the effects of varying 
these choices. 

First, we consider the effects of using coordinate time t 
in place of r. In Fig. ED we plot the natural loga- 




T(M) 

FIG. 16: Natural logarithms of the ellipsoid axes vs. r for 
rescaling time T = 100 M (dark) and time T = 50 M (light). 

rithm of the largest ellipsoid axis log[ri(r)] vs. r together 
with log[ri(f)] vs. t. (We use the unrescaled deviation 
vector approach for simplicity, since the Jacobian ap- 
proach requires a new Jacobian matrix for each coordi- 
nate change.) The exponents are A r = 5.05 x 10~ 3 M -1 
and A* = 2.51 x 10~ 3 M^ 1 , implying Lyapunov timescales 
of t \ = 1.98 x 10 2 M and t x = 3.98 x 10 2 M . The 
average value of dt/dr over the orbit is 2.06, whereas 
t\/r\ = 2.01, so the relationship 

Hi) 

discussed in the introduction is well satisfied. 

Second, we discuss the effects of varying the rescal- 
ing time T. We find that choosing T to be a moderate 
fraction of the shortest Lyapunov timescale (correspond- 
ing to the largest Lyapunov exponent) works best, giving 
each axis enough time to grow before rescaling while still 
keeping the negative exponents from underflowing and 
preventing the largest axis from dominating. Rescaling 
times between 50 M and 100 M work best for the sys- 
tems we consider, which have Lyapunov timescales rang- 
ing from 10 2 M to 10 3 M. A comparison of results for 
T = 50 M and T = 100 M appears in Fig. [H 

Third, we compare the projected norm used here to a 
naive Euclidean norm for determining the length of the 
phase-space tangent vectors £i. As shown in Fig. El 
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FIG. 17: The natural logarithm of the largest ellipsoid axis 
vs. r for the Euclidean norm (top) and the projected norm 
from Sec. HVaI (bottom) . 



even using a 12-dimensional Euclidean norm changes the 
resulting exponent very little (approximately 15% in this 
example). Given its conceptual advantages, we choose 
to use the projected norm with the confidence that the 
Lyapunov exponent order of magnitude is robust. 



VIII. CONCLUSIONS 

A spinning test particle, as described by the Papa- 
petrou equations, appears to be chaotic in Kerr space- 
time, with maximum e-folding timescales of a few x 
10 2 M . The applicability of this result is limited by three 
main factors: (1) chaos appears only for physically unre- 
alistic values of the spin parameter; (2) other effects, such 
as tidal coupling, may be important for some astrophysi- 
cal systems, violating the pole-dipole approximation im- 
plicit in the Papapetrou equations; and (3) we neglect 
gravitational radiation. The third limitation is not fatal, 
since the radiation timescales can be long enough that 
chaos, if present in the conservative limit, would have 
time to manifest itself in the gravitational radiation of 
extreme mass-ratio systems. 

In the unphysical S = 1 limit, the Lyapunov expo- 
nents exhibit characteristics expected of a Hamiltonian 
system, appearing in pairs ±A (Sec. IVI Bp . There are 
zero Lyapunov exponents which correspond to the con- 
stants of the motion, but the other exponents are in gen- 
eral nonzero. (For the Kerr orbits considered in this 
paper, we find that two of the three independent ex- 
ponents are nonzero, as illustrated in Fig. 0) Typical 
orders of magnitude for the largest Lyapunov exponents 
are a few x 10 _3 M _1 for unphysical spins (S = 1). For 
physically realistic spin parameters (Sec. IVlfCl) . we find 
that A max < few x 10 _7 M _1 , corresponding to e-folding 
timescales of a few x 10 6 M. Even this bound appears to 
be limited only by the total integration time; in all phys- 
ically realistic cases considered, A max is indistinguishable 
from zero (using S = integrations as a baseline) . 

From the perspective of gravitational radiation detec- 



tion, our most important conclusion is that chaos seems 
to disappear for physically realistic values of S, i.e., val- 
ues of S for which the test-particle approximation and 
hence the Papapetrou equations are valid. We are unable 
to comment on the dynamics of comparable mass-ratio 
binaries, since such systems are not accurately modeled 
by the Papapetrou equations, but for extreme mass-ratio 
binaries it appears unlikely that chaos will present a prob- 
lem for the calculation of theoretical templates for use in 
matched filters. A more thorough exploration of param- 
eter space is needed to reach a firmer conclusion [? ] . 
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APPENDIX: FULL JACOBIAN 

For reference, we list the derivatives needed to calcu- 
late the full Jacobian matrix. 

From Sec. IV Bl we have the following: 

— = N[ Pag ^ + w% 
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Now we simply apply the product rule many times: 
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